Thanks @Mashin6, that looks great. One immediate thought is that the sub-divided parts of a bbox really need to have a small degree of overlap to ensure that all objects that cross the boundary lines are appropriately captured. I like the way you've coded the bbox division - would it be possible to extend that elegant simplicity to include an additional parameter to control proportional overlap across boundaries? I usually use a default overlap of 5%, and once had code somewhere that showed that it did make a difference. (That was, however, from an older version of overpass, so maybe no longer valid?)

Yeah, I guess one can add a fudge factor that would extend the newly made boundaries by a bit.

split_bbox <- function(bbox, grid = 2, ff = 0.05) {
    ...
    b <- matrix(c(xmin + ((i-1) * dx),
                  ymin + ((j-1) * dy),
                  xmin + (i * dx) + (dx * ff),
                  ymin + (j * dy) + (dx * ff)),
                nrow = 2,
                dimnames = dimnames(bbox))
    ...
}

But we should not extend past the original bbox, so there must be a switch for the last edge. e.g.

split_bbox <- function(bbox, grid = 2, ff = 0.05) {
    ...
    b <- matrix(c(xmin + ((i-1) * dx),
                  ymin + ((j-1) * dy),
                  xmin + (i * dx) + (dx * ff * (if (i == grid) 0 else 1)),
                  ymin + (j * dy) + (dx * ff * (if (j == grid) 0 else 1))),
                nrow = 2,
                dimnames = dimnames(bbox))
    ...
}

Though I don't know how much it matters. Overpass seem to return also objects that are directly positioned on the bbox edge. For example this node has coordinates 41.3084017, -72.9262955 and is returned by this query node(41.3084017, -72.9262955, 41.3085, -72.925);. OSM works with 7 decimal places and coordinates with more precision get rounded. So potential problem could be when the new bbox after splitting has coordinates with more than 7 decimal places and it gets rounded in the wrong direction. Then at least in theory it should be sufficient to always round up the 7th digit of the x and y max coordinates like:

    b <- matrix(c(xmin + ((i-1) * dx),
                  ymin + ((j-1) * dy),
                  ceiling( (xmin + (i * dx)) * 10^7) / 10^7,
                  ceiling( (ymin + (j * dy)) * 10^7) / 10^7),
                nrow = 2,
                dimnames = dimnames(bbox))

Great, thanks. The border thing I've previously addressed with a final line to trim all sub-boxes to lie within the original box, so that only internal boundaries get fudged. And the previous observations, which I can unfortunately no longer find, were way/rel objects which were excluded because they extended beyond the bbox. Those should be included by overpass recurse up, but i definitely had some cases where that didn't happen. But maybe overpass updates since then would ensure that this would no longer happen, so maybe it's not even necessary? I'll play around myself a bit, and maybe end up asking you to undo what i just asked you to do - but even if so, it will have been worthwhile, because it is important to explicitly note that recurse up (hopefully) ensures that all objects with any node within a bbox are included in a query. I'll get back to you asap. Thanks for the work!